This document shows how to calculate candidate Pillars and their constituent Indicators at Census Block Group, Tract, and Service Area/ Municipality level. Here, we load the data we pre-processed in the previous exercise
b <- read_sf("../data/boundary.geojson")
bg.l <- read_sf("../data/blockgroups_long.geojson")
bg.w <- read_sf("../data/blockgroups_wide.geojson")
tr.l <- read_sf("../data/tracts_long.geojson")
tr.w <- read_sf("../data/tracts_wide.geojson")
pl.w <- read_sf("../data/place_wide.geojson")
pl.l <- read_sf("../data/place_long.geojson")
lead <- read_sf("../data/parcels_lead_service_lines.geojson")
lead.c <- sf::st_centroid(lead)
calculate_volume <- function(hh_size = 4, gpcd = 50){
vol = hh_size * gpcd * 30
return(vol)
}
calculate_bill <- function(volume = 6000){
fixed_charges = 7.63 + 9.85 + 1.8
vol_rate = 5.76 + 2.71
bill = fixed_charges+ (vol_rate*volume)/748.052
return(bill)
}
gradeAffordability <-function(HBI,PPI){
grade <- NA
grade <- case_when(
HBI < 7 & PPI < 20 ~ "Low Burden",
HBI < 7 & PPI>=20 & PPI <35 ~ "Moderate-Low Burden",
HBI < 7 & PPI >=35 ~ "Moderate-High Burden",
HBI >= 7 & HBI < 10 & PPI >= 35 ~ "High Burden",
HBI >= 7 & HBI < 10 & PPI >= 20 & PPI < 35 ~ "Moderate-High Burden",
HBI >= 7 & HBI < 10 & PPI < 20 ~ "Moderate-Low Burden",
HBI >= 10 & PPI >= 35 ~ "Very High Burden",
HBI >= 10 & PPI >= 20 & PPI < 35 ~ "High Burden",
HBI >= 10 & PPI < 20 ~ "Moderate-High Burden"
)
return(grade)
}
save.image("../cache/components_prep.RData")
load("../cache/components_prep.RData")
We adopt the approach in the figure below, as elaborated in this report.
HBI/PPI Matrix
#HBI Calculation
# Place
pl.w <- pl.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
pl.w <- pl.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
pl.w <- pl.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
pl.w <- pl.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
pl.w <- pl.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
pl.w <- pl.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
pl.w <- pl.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
pl.w <- pl.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))
pl.w$HBI_size_avg <- 100*pl.w$bill_size_avg/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_1 <- 100*pl.w$bill_size_1/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_2 <- 100*pl.w$bill_size_2/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_3 <- 100*pl.w$bill_size_3/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_4 <- 100*pl.w$bill_size_4/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_5 <- 100*pl.w$bill_size_5/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_6 <- 100*pl.w$bill_size_6/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$HBI_size_7 <- 100*pl.w$bill_size_7/(pl.w$hh_income_upper_limit_Q1E/12)
pl.w$PPI <- 100*(pl.w$inc_poverty_countE-pl.w$inc_poverty_count_gr2E)/pl.w$inc_poverty_countE
pl.w$aff_grade_avg <- gradeAffordability(pl.w$HBI_size_avg,pl.w$PPI)
pl.w$aff_grade_1 <- gradeAffordability(pl.w$HBI_size_1,pl.w$PPI)
pl.w$aff_grade_2 <- gradeAffordability(pl.w$HBI_size_2,pl.w$PPI)
pl.w$aff_grade_3 <- gradeAffordability(pl.w$HBI_size_3,pl.w$PPI)
pl.w$aff_grade_4 <- gradeAffordability(pl.w$HBI_size_4,pl.w$PPI)
pl.w$aff_grade_5 <- gradeAffordability(pl.w$HBI_size_5,pl.w$PPI)
pl.w$aff_grade_6 <- gradeAffordability(pl.w$HBI_size_6,pl.w$PPI)
pl.w$aff_grade_7 <- gradeAffordability(pl.w$HBI_size_7,pl.w$PPI)
pl.bill<- pl.w %>% st_drop_geometry() %>% select(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7) %>%
pivot_longer(cols = c(bill_size_avg,bill_size_1,bill_size_2,bill_size_3,bill_size_4,bill_size_5,bill_size_6,bill_size_7),names_to="hh_size",names_prefix="bill_size_",values_to="bill")
pl.HBI<- pl.w %>% st_drop_geometry() %>% select(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7) %>%
pivot_longer(cols = c(HBI_size_avg,HBI_size_1,HBI_size_2,HBI_size_3,HBI_size_4,HBI_size_5,HBI_size_6,HBI_size_7),names_to="hh_size",names_prefix="HBI_size_",values_to="HBI")
pl.aff<- pl.w %>% st_drop_geometry() %>% select(PPI,hh_income_upper_limit_Q1E,aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7) %>%
pivot_longer(cols = c(aff_grade_avg,aff_grade_1,aff_grade_2,aff_grade_3,aff_grade_4,aff_grade_5,aff_grade_6,aff_grade_7),names_to="hh_size",names_prefix="aff_grade_",values_to="aff_grade")
pl.aff <- pl.aff %>% left_join(pl.HBI,by="hh_size") %>% left_join(pl.bill,by="hh_size") %>% select(hh_size,bill,aff_grade,HBI,hh_income_upper_limit_Q1E,PPI)
print(paste0("The Average Household Size for Naperville is ",
round(pl.w$hh_size_avg_overallE,2)))
[1] “The Average Household Size for Naperville is 2.79”
print(paste0("At 50 GPCD, this corresponds to a monthly water volume of ",
calculate_volume(hh_size=pl.w$hh_size_avg_overallE)))
[1] “At 50 GPCD, this corresponds to a monthly water volume of 4185”
print(paste0("This corresponds to an water and sewer bill of $",
round(pl.w$bill_size_avg,2)))
[1] “This corresponds to an water and sewer bill of $66.67”
print(paste0("The Upper limit of the lowest quintile of annual household income for Naperville is $",
pl.w$hh_income_upper_limit_Q1E,
" or $",
round(pl.w$hh_income_upper_limit_Q1E/12,2),
" monthly"))
[1] “The Upper limit of the lowest quintile of annual household income for Naperville is $57958 or $4829.83 monthly”
print(paste0("The Poverty Prevalence Indicator (PPI) for Naperville overall is % ",
round(pl.w$PPI,1)))
[1] “The Poverty Prevalence Indicator (PPI) for Naperville overall is % 9.4”
print(paste0("By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7"))
[1] “By the HBI/PPI Matrix, the Water and Sewer Cost Burden in Naperville overall is graded below, assuming an average housheold size as well as household sizes 1-7”
datatable(pl.aff) %>%
formatCurrency(c('bill','hh_income_upper_limit_Q1E'),'$') %>%
formatRound(c('HBI','PPI'),2)
Below we visualize HBI, PPI, and the overall affordability grade as calcualted at the tract level.
# Tract
# Place
tr.w <- tr.w %>% mutate(bill_size_avg = calculate_bill(volume = calculate_volume(hh_size=hh_size_avg_overallE)))
tr.w <- tr.w %>% mutate(bill_size_1 = calculate_bill(volume = calculate_volume(hh_size=1)))
tr.w <- tr.w %>% mutate(bill_size_2 = calculate_bill(volume = calculate_volume(hh_size=2)))
tr.w <- tr.w %>% mutate(bill_size_3 = calculate_bill(volume = calculate_volume(hh_size=3)))
tr.w <- tr.w %>% mutate(bill_size_4 = calculate_bill(volume = calculate_volume(hh_size=4)))
tr.w <- tr.w %>% mutate(bill_size_5 = calculate_bill(volume = calculate_volume(hh_size=5)))
tr.w <- tr.w %>% mutate(bill_size_6 = calculate_bill(volume = calculate_volume(hh_size=6)))
tr.w <- tr.w %>% mutate(bill_size_7 = calculate_bill(volume = calculate_volume(hh_size=7)))
tr.w$HBI_size_avg <- 100*tr.w$bill_size_avg/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_1 <- 100*tr.w$bill_size_1/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_2 <- 100*tr.w$bill_size_2/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_3 <- 100*tr.w$bill_size_3/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_4 <- 100*tr.w$bill_size_4/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_5 <- 100*tr.w$bill_size_5/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_6 <- 100*tr.w$bill_size_6/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$HBI_size_7 <- 100*tr.w$bill_size_7/(tr.w$hh_income_upper_limit_Q1E/12)
tr.w$PPI <- 100*(tr.w$inc_poverty_countE-tr.w$inc_poverty_count_gr2E)/tr.w$inc_poverty_countE
tr.w$aff_grade_avg <- gradeAffordability(tr.w$HBI_size_avg,tr.w$PPI)
tr.w$aff_grade_1 <- gradeAffordability(tr.w$HBI_size_1,tr.w$PPI)
tr.w$aff_grade_2 <- gradeAffordability(tr.w$HBI_size_2,tr.w$PPI)
tr.w$aff_grade_3 <- gradeAffordability(tr.w$HBI_size_3,tr.w$PPI)
tr.w$aff_grade_4 <- gradeAffordability(tr.w$HBI_size_4,tr.w$PPI)
tr.w$aff_grade_5 <- gradeAffordability(tr.w$HBI_size_5,tr.w$PPI)
tr.w$aff_grade_6 <- gradeAffordability(tr.w$HBI_size_6,tr.w$PPI)
tr.w$aff_grade_7 <- gradeAffordability(tr.w$HBI_size_7,tr.w$PPI)
hbi <- select(tr.w,
HBI_size_avg,
HBI_size_1,
HBI_size_2,
HBI_size_3,
HBI_size_4,
HBI_size_5,
HBI_size_6,
HBI_size_7,
)
ppi <- select(tr.w,
PPI)
aff <- select(tr.w,
aff_grade_avg,
aff_grade_1,
aff_grade_2,
aff_grade_3,
aff_grade_4,
aff_grade_5,
aff_grade_6,
aff_grade_7)
m1 <- mapview(aff,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m2 <- mapview(hbi,burst=TRUE,col.regions=viridis::magma(7)) + mapview(b,alpha.regions=0,col.regions="green",stroke=TRUE,lwd=3,color="green",layer.name="Municipal Boundary")
m3 <- mapview(tr.w,zcol="hh_income_upper_limit_Q1E",col.regions=viridis::plasma(7),layer.name="Lower Quintile Income") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m4 <- mapview(ppi,burst=TRUE) + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
sync(m1,m2,m3,m4)
Below, we simulate cutoff rates as a probabilistic function of HBI and PPI, and Customer Assistance participation rates as a probabilistic function of cutoff rates and white population share. We simulate participation in conservation/ efficiency appliance incentive programs as a function of homeownership rates and rate of population with at least a bachelors degree.
In reality, one would hope to use actually utility records of participation in these programs at the address level, or aggregated to census tracts.
tr.w$white_perc <- tr.w$pop_race_whiteE/tr.w$pop_race_countE
P1 <- select(tr.w, NAME, GEOID, aff_grade_avg, HBI_size_avg, PPI, hh_income_upper_limit_Q1E, hh_size_avg_overallE,white_perc, pop_owner_occupied_unitsE, pop_rental_unitsE,
pop_educ_attainment_countE, pop_educ_bachelorE, pop_educ_masterE, pop_educ_docE, pop_educ_profE)
set.seed(20212002)
#Cutoffs and CAP
z <- 5 -0.03* P1$HBI_size_avg - 0.02*P1$PPI - 0.04 * P1$HBI_size_avg*P1$PPI
P1$Cutoff_Perc = 100*(1/(1+exp(z)))
P1$CAP_Perc <- P1$Cutoff_Perc * (1-0.3*(1-P1$white_perc))
#Conservation (ever)
P1$bachAbove <- 100*(P1$pop_educ_bachelorE + P1$pop_educ_masterE + P1$pop_educ_docE + P1$pop_educ_profE)/P1$pop_educ_attainment_countE
P1$own_rate <- 100*P1$pop_owner_occupied_unitsE/(P1$pop_rental_unitsE+P1$pop_owner_occupied_unitsE)
z <- 4 - 0.01 * P1$bachAbove - 0.02 * P1$own_rate
P1$Incentive_rate <- 100*(1/(1+exp(z)))
m5 <- mapview(P1,zcol="Cutoff_Perc", layer.name="Account cutoff (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m6 <- mapview(P1,zcol="CAP_Perc", layer.name="CAP participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
m7 <- mapview(P1,zcol="Incentive_rate", layer.name="Efficiency Incentive participation (%)") + mapview(b,alpha.regions=0,col.regions="red",stroke=TRUE,lwd=3,color="red",layer.name="Municipal Boundary")
sync(m5,m6,m7)
We model all water quality events as basically random occurrences with probabilities between 0 and 10% and coverage of 50-100% at the tract level
set.seed(32896)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$Interrup <- event*coverage
set.seed(328435396)
prob <- runif(length(P1$GEOID), min=0,max=0.1)
coverage <- runif(length(P1$GEOID),min=50, max=100)
event <- rbinom(length(P1$GEOID),1, prob=prob)
P1$boil_rate <- event*coverage
This would be a Utility-Level measure only, and is available from EPA SDWIS
In this case, the API call for Naperville involves the PWSID IL0434670:
https://enviro.edap-cluster.com/cache/query/enviro3/efservice/VIOLATION/PWSID/IL0434670/IS_HEALTH_BASED_IND/Y/ROWS/0:100/JSON
In this case, there are no health-based violations in the past 10 years, so the value of the indicator is 0.
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.